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INTRODUCTION 

There  exist  a  large  number  of  planetary  boundary  layer  models  each  designed 
for  specific  purposes,  (Deardorff, .1974;  Mahrt  and  Lenschow,  1976;  Stull,  1976; 
Heidt,  1977;  Yamada,  1979;  Tennekes  and  Driedonks,  1981;  Chen  and  Cotton, 

1983) .  Many  of  these  consider  the  top  of  the  boundary  layer  a  material 
surface.  Some  consider  the  top  of  the  boundary  layer  to  be  coincident  with 
the  inversion,  and  consider  entrainment  across  its  interface.  Some  are  multi  - 
level,  some  are  bulk,  single  -  level  models.  The  various  models  are  of  one, 
two  and  three  dimensions.  The  greater  the  number  of  dimensions,  the  greater 
the  computational  complexity  . 

It  is  possible  to  reduce  the  computational  complexity,  and  still  not 
eliminate  the  three  -  dimensionality  completely  by  using  a  variation  of  the 
"momentum  integral"  method  (Schlichting,  1968).  By  means  of  this  approach,  the 
vertical  structure  of  several  of  the  variables  is  specified  and  incorporated 
into  the  vertically  integrated  equations.  In  this  way,  the  model  becomes  two  - 
dimensional  in  the  horizontal,  and  vertical  variations  are  incorporated  into 
various  coefficients.  The  method  was  introduced  into  meteorology  by  Charney 
and  Eliassen  (1979),  who  derived  the  highly  successful  "equivalent  barotropic 
model".  The  method  has  also  been  used  in  studying  nocturnal  drainage  flow 
(Manins  and  Sawford,  1979). 

In  this  study,  one  of  the  guiding  principles  was  that  computers  of  limited 
power  would  be  available  to  us.  Therefore  we  decided  at  the  outset  to  attempt 
to  model  the  atmospheric  circulation  in  the  desert  planetary  boundary  layer  by 
means  of  a  vertically  -  integrated,  parameterized  model. 
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We  consider  a  model  of  the  planetary  boundary  layer  (depth  approximately 
lkm) .  This  layer  is  itself  divided  into  a  surface  layer  (20m) ,  a  transition 
layer,  which  at  certain  times  and  places  becomes  very  well  mixed  (and  is 
frequently  called  the  "mixed  layer")  and  an  inversion  layer.  (See  Figure  1) . 

We  shall  not  assume  that  the  transition  layer  is  always  well  mixed. 

Very  often,  the  top  of  the  planetary  boundary  layer  is  capped  by  an 
inversion.  This  is  particularly  true  in  many  desert  and  semi  -  desert  regions. 
In  Israel  for  example,  there  are  222  days  per  year,  on  the  average,  of  mid-day 
inversions  just  about  100  km  north  of  the  beginning  of  the  Negev  desert.  (Shaia 
and  Jaffe  1976)  As  we  go  farther  south,  closer  to  the  center  of  the  Hadley 
cell,  the  frequency  of  occurrence  of  inversions  is  probably  even  greater.  When 
the  capping  inversion  exists,  and  when  convection  occurs  below  the  inversion, 
the  inversion  changes  height  due  to  upward  and  corresponding  downward  fluxes 
through  it  by  turbulence.  The  inversion  height  changes  affect  the  dust 
concentration.  These  processes  have  to  be  modelled.  Further,  the  processes 
which  we  wish  to  model  are  on  such  a  scale  that  fairly  high  resolution  is 
needed  on  the  order  of  10-20  km  in  the  horizontal,  over  a  region  approximately 
300x600  km.  If  then  the  vertical  resolution  above  the  surface  layer  is  to  be 
very  fine  -  say  100m  -  the  computation  time  for  solution  of  only  the  boundary 
layer  mesoscale  equations  may  become  prohibitive.  Thus,  in  spite  of  the  fact 
that  the  optimum  mesoscale  model  must  be  three-dimensional,  we  expect  to  gain 
valuable  insights  with  a  vertically  parameterized  two-dimensional  nonsteady 


We  shall  concern  ourselves  with  a  form  of  the  primitive  equations  derived 


by  ensemble  averaging  over  a  horizontal  area  4x&y  which  is  large  enough  to 
contain  the  sub-grid  scale  phenomena,  but  small  enough  to  be  a  fraction  of  the 
mesoscale  system.  We  define 

At  A ^ 

O'  —  ~C*  A- 

whereo^is  any  scalar  variable. 

With  the  above  definitions,  the  appropriate  equations  are,  approximately 
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In  the  above  set  of  equations,  the  unbarred  variables 
u,v,w,Ug,Vg,e,F,P,q,c,  are  all  mean  values  according  to  Equation  (1). 

We  have  used  the  Boussinesq  assumption,  and  the  variables  are  defined  as 
follows: 


u,v,w 


u'  ,v'  ,w' 

U  .V 

9  9 
f 

e 

P 

F 

q 

c 

e* ,q' »c* 
x,y,z,t 


components  of  the  wind  vector 

turbulent  components  of  the  wind  vector 

geostrophic  wind  components 

Coriolis  parameter 

potential  temperature 

air  density 

net  radiation  flux 

specific  humidity 

dust  concentration 

turbulent  fluctuations  of  0,q,c 

space  variables  and  time 


Cp  =  specific  heat  of  air  at  constant  pressure 

SL(r)  =  fall  velocity  of  particle  of  radius  r 

Equations  (2)  and  (3)  are  the  horizontal  equations  of  motion.  Equation  (4)  is 
the  continuity  equation.  Equation  (5)  is  the  thermodynamic  energy  equation. 
Equation  (6)  is  the  moisture  equation,  Equation  (7)  is  the  dust  concentration 
equation. 

The  lower  boundary  condition  is 

w  =  wT  =V-v  zT  at  z  =  zT  =  terrain  height  (8) 

Here  "V  =  t  u  +  jj  v  =  horizontal  wind  vector,  l  and  are  unit  vectors  in  the  E 
and  N  directions  respectively  and  V  is  the  gradient  operator  in  2  -  dimensions 

V  =  l— +  A\  . 

In  the  present  investigation,  we  assume  that  any  condensed  moisture  stays 
in  the  air.  Thus  we  do  not  treat  clouds  or  precipitation  explicitly  in  this 
model  (but  implicit  predictions  are  possible) . 

In  order  to  expedite  deriving  appropriate  expressions  for  predicting 
inversion  height,  we  first  derive  inversion  "interface"  conditions.  The 
results  are  essentially  the  upper  boundary  conditions  .  Let  h(x,y,t)  be  the 
inversion  height.  Let  i  be  a  small  layer  of  constant  thickness  above  the 


inversion  level.  In  Mahrt  and  Lenschow  (1976)  this  is  called  the  turbulent 
inversion  layer,  which  is  sufficiently  thin  so  that  terms  of  0(£)  may  be 


neglected  in  the  integrated  equations. 


Define 


L<Y)  =  -L-  \  ( o,\  dz 
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(9) 


Let 


4  O'  =  +  J 

=  jump  in  of  across  inversion 


(10) 


By  Leibniz's  rule, 
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ia.  As r  =  fa. 


where 
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>jy 


(11) 


We  apply  the  operator  Equation  (9)  to  Equations  (2) ,  (3) ,  (5) ,  (6) , 
allow  to  approach  zero,  and  obtain. 


~  AU  +  +  ^CuvyzL-  -  -  (u'  */') 
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Ok  -  44  4  4  C-V*)^  =  -  (  V  W') 


(12) 


(13) 
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(14) 


In  the  above  derivations,  we  have  assumed  that 


=  (ii4  w-\ 

Each  of  the  above  equations  can  be  viewed  as  equations  for  the  vertical 


eddy  fluxes,  or  as  prediction  equations  for  h  if  these  eddy  fluxes  are  known. 
The  quantity  (  -  w^) ,  which  represents  the  motion  of  the  air 

relative  to  that  of  the  inversion,  is  called  the  "entrainment  velocity",  wg 
is  the  larger  scale  vertical  velocity  at  the  interface.  The  terms 
involving  “ih/gx  and  ^h/gy  are  usually  omitted  in  derivations  of  these  interface 
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conditions,  since  most  inversion  height  models  assune  horizontal  homogeneity. 


Our  approach  will  be  to  integrate  the  system  of  Equations  (2)- (7) 
inclusive,  with  respect  to  z  from  z  *  k  =  constant  =  height  of  surface  layer, 
to  z  =  h  =  inversion  height,  then  to  introduce  modelling  assumptions  for  all 
the  variables,  i.,  e.,  to  specify  their  variations^with  height.  If  we  do  this, 
it  becomes  possible  to  express  all  of  the  jump  quantities  in  Equations 
( 12 ) — ( 16 )  inclusive  in  terms  of  their  values  at  specific  levels. 


We  assume  that  the  surface  layer  is  neutral,  so  that  we  may  use  a  constant 
flux  profile 

Vfr^fc.fc)  <17> 

•ft, 

where  ^  *  von  Karman  constant  =  0.4,V*  =  friction  velocity, 

2  2  1/2 

zQ  =  roughness  parameter.  Here  V*  =  (u*  +v*  )  , 


where  u*and  v*  are  defined  below. 

Transition  Layer: 

We  assume 

u  (x,y,z,t)  =  A(z)  u  (x,y,t).  *  r 
v  (x,y,z,t)  »  B(z)  v  (x,y,t) 

where 

(£}  -  77— r  C  ^  ^  1 

[  ^  **  'x)  J  -ft. 

At  the  level  z  *  k,  the  two  expressions  (17)  ,  (18)  must  match 

U  (4)  =  H*_  A.  (A±ii)  =  a^C) 


Therefore 


potential  Temperature 


Surface  Layer 


Assuming  that  the  turbulent  flux  of  heat  is  also  constant  with  height  in  the 
surface  layer,  we  find 


e  (x,y,z,t)  =  eGR  ♦  Mk-  eGRi 


A  fir) 


We  derive  a  modelling  approximation  for  potential  temperature  in  regions 
where  a  nighttime  inversion  exists,  and  where  surface  heating  during  the  day 
leads  to  convection  and  turbulent  mixing.  In  Figure  1,  an  inversion  exists  in 
the  early  morning.  The  potential  temperature  at  z  =  k  is  *  6(x,y,k,0) . 

The  potential  temperature  increases  linearly  with  height  according  to 


0(x,y,z,t)  =  0kI  +  Y(0) (z-k) 


up  to  z=h,  and  linearly  from  there  up  to  z  =  ht$,  with  a  lapse  rate  1^0)  within 
the  inversion  layer.  Here  f(0)  is  also  the  lapse  rate  above  the  inversion 


©(x,y,hfi,t)  =  =  ©kI  +  r«))(z-k)  +  Vinv(0)f 


(23) 


It  is  assumed  that  heating  destabilizes  the  lapse  rates,  so  that,  at  some  later 
time 

0(x,y,z,t)  =  ©k  +  Wt)  (z-k) 

6(x,y,h,t)  =  ©h  =  ©R  +  f(t)  (h-k) 

e(x,y,h+^,t)  =  ©h  +i  =  ek  +  Ht)  (h-k)  +  rinv(t)£ 

It  is  assuned  that  Equations  (23)  and  (25)  will  ] 

.  In  that  case, 

4  e  =  (BkI-ek)  +  tr(0)-r(t)i  (h-k)  +  rinv(0)  S  (26) 

From  Equation  (26),  we  see  that 

and,  since 

from  first  of  Equations  (25),  we  have 


(24) 


(25) 


field  the  same  result  for 


(27) 


This  is  similar  to  the  result  obtained  by  Tennek.es  and  Driedonks  (1981)  for 


a  well-mixed,  horizontally  homogeneous  layer,  i.e.,  "the  magnitude  of  the 
temperature  junp  0  increases  in  two  ways:  it  increases  as  the  height  of  the 
mixed  layer  increases,  and  decreases  if  the  boundary  layer  warms  up".  Here  we 
have  not  assumed  well-mixedness  (©h  #©m)  or  horizontal  homogeneity,  so 
that  Equation  (27)  is  more  general  than  was  realized. 

In  the  case  of  flow  over  water,  the  lapse  rate  'f(t)  is  replaced  by  the 
appropriate  expression  for  a  marine  environment,  say  Y^(t). 


Moisture 


We  follow  the  same  formulation  as  for  potential  temperature. 


q(x,y,z,t)  =  qk(x,y,t)  +  $(t) (z-k) 


(28) 


q(x,y,z,0)  =  qRI  +5(0)  (z-k) 


(29) 


4q  =  (qkI  -  qk)  +  [*(0)  -  J(t)]  (h-k)  +  Sinv(0)S 


(30) 


In  the  case  of  flow  over  water,  the  lapse  rate  of  moisture  3(t)  is  replaced 
by  its  appropriate  value  over  water,  say  $^(t). 

Dust  p 

We  assume  that  the  dust  concentration  is  given  by 

C(x,y,z,t)  =  D(z)  Ck(x,y,t)  (31) 

So  that 

4  c  =  [D(hf^)  -  D(h)  1  Ck(x,y,t)  (32) 

If  we  make  the  further  assunption  that  there  is  no  dust  at  the  top  of  the 
inversion  layer,  i.e.,  at  z  =  h  +£  ,  then  D(hff)  =  0,  (See  Carlson  and 
Prospero,  1972,  concerning  top  of  Saharan  dust  layer),  and 


4  c  =  -  D  (h)  C.  (x,y,t) 


(33) 


Surface  Parameterizations 


For  the  momentum,  we  use 


(u’w*)k  =  -Cd  (A(k)l  |*|  u 


(v'w')k  =  -Cd  [  B(k)  3  I*!  v 


where  C ,  is  the  turbulent  transfer  coefficient  for  momentum 
a  * 

c .  =  5xl0'4  V1^2  (V  in  ms'1)  (Berkofsky,  1982) 

a 

For  the  convective  heating, 


(",S,)K  =CH0  A(K)  U  (9GR  *  9k> 


where  0u,R  is  the  potential  temperature  at  the  ground,  is  the 


For  the  moisture, 


(Wq')k  -  A(k)  u  [  qsat  (e-GR)  -  qkI 
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where  q__.  (©  )  means  saturation  mixing  ratio  at  temperature 

s<3  z  or 


0__,  W_n  =  ground  soil  moisture,  W.  =  potential  saturation 

CjH  CjK  k 


value  of  W. 


For  the  dust. 


(w'c')k  =  Cd  A(k)  u  (C™  ~  C ^ 


GR  k' 


(37) 


where  CGR  is  the  dust  concentration  in  a  thin  layer  near  the  ground. 


Ground  Albedo 


or  = 

gr 


a  +  b  ^  ,b<0 

w. 


(38) 


a,b  constants 


In  the  above,  WGR,  ©GR,  CGR  must  be  predicted. 


Radiation 

The  formulation  of  short  and  longwave  energy  exchange  in  the  air  and  at  the 
ground  is  based  on  Gates  et  al  (1971).  The  system  is  applied  at  three  levels: 
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ground,  top  of  surface  layer,  and  top  of  transition  layer:  Thus  the  temperatures 
at  these  levels  (z  =  0,  z=k,  z  =  h)  are  needed,  as  is  the  temperature  at  2m. 
All  of  these  are  obtainable  within  the  framework  of  the  model  by  converting 


potential  temperatures  as  given  to  .temperatures  using  T  =  8(  -*-)* 

P. 


Longwave  radiation  balance  (  positive  if  upwards)  at  each  of  the  three 
levels  is  defined  by  F^,  and  RGR  as  follows: 

*  o->U  \fT% T (u*  - +  f (T^-Tp  Li |  +  O. tearful) 


^Ta  [o-t  h(u^)  -  O-l]  4  ^ 


where 


i  + 


C*  =  -XT) 


T  =  temperature  at  2m. 

3 


<r  =  Stef an- Boltzmann  constant 


(39) 


(40) 


u.  is  the  effective  water  vapor  content  between  ground  level  and 


i,  and  is  obtained  from 


(41) 


where  p^r  pressure  at  the  ground  p^  =  pressure  at  level  i. 


i '  *■  (±)' , 


(42) 


where  is  the  mixing  ratio  at  the  ground  and  \  =  constant  =  2.92 
(Smith,  1966) . 


For  levels  k,  h,o*»,  we  obtain 


X+X  X+i 

h*  -  P* 


*  +» 


Up, 


-  4 


6-8 


*c 


-  4 


»» 


1r  &*■ 


? 

Sr. 


-IP, 


x+* 

* 


x+r 


x*> 


*  +* 
-  h 


<50 


>v+> 


Note:  if  we  assume  u^*  =  0,'T  (uk*)=  1  ^ (u*-uk*)  =T(iO. 

4  4  4 

Also,  if  we  assume  Tk  »  (Tfl  -Tk  ) ,  we  obtain 
a  simple  form  for  Rk>  We  have  not  made  use  of  these. 


(43) 


The  effect  of  CC>2  absorption  is  taken  into  account  by  the 
coefficients  0.736  and  the  term  O.o'j  *\  (u^)  -  0.1.  The  former  value 
actually  applies  at  600  mb  in  the  Mintz  -  Arakawa  model,  but  we  use  this 
for  h  and  k,  which  are 
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very  much  lower  than  600  mb.  The  effect  of  the  error  is  not  known. 


Shortwave  radiation  (positive  if  downwards)  is  given  by 

^  r  S* 1 1-  [fwi  -  M.J  )jet  ’”1 

%  •  <  -<  {  I-  )S»c  j]*'i”  ] 

<  .S 

*>6-*  ~  +  5* ^ 

$4*  =  ^  fl- (VO 

(\-  <Yo<V*<0 

where  SQa  =  part  of  the  solar  radiation  subject  to  absorption  = 
SQ  x  0.651  cos  Z 

SQs  =  part  of  the  solar  radiation  subject  to  scattering  = 

S0  x  0.„h9  cos  Z 


©Tgr  =  ground  albedo, 

Cy-0  =  min  ^1,  0.085  -  0.247  log  10  [  —  cosZ 

=  sky  albedo,  ( 

Z  =  zenith  angle  of  the  sun 

Sq  =  solar  constant 


The  radiation  balance  =  outgoing  -  incoming,  i.e.  negative 
radiation  balance  at  the  ground  means  input  of  radiative  heat  to  the 


surface. 


Closures 


At  this  point,  we  define  several  closure  formulas,  and  will  indicate 
later  where  and  why  they  are  used. 


(ur’e')*  =-A. 

\  =  fl*  (*)  (  £1  ' 


Cov»st*VNk  ,  \  ^ 


We  are  now  in  a  position  to  derive  the  final  form  of  the  equations. 

The  procedure  is  to  apply  the  averaging  Equation  (19)  to  Equations  (2) — (V) 

inclusive  and  use  the  interface  conditions  Equations  (12)-(16)  inclusive  to 

eliminate  the  fluxes  at  inversion  height,  thus  obtaining  the  transition 

layer  equations.  If  we  do  this  in  all  of  the  equations,  we  are  left 

"bJi. 

without  a  prediction  equation  for  .  However  if  we  use  the  first  of 
Equations  (44)  ,  which  states  that  the  turbulent  flux  of  heat  at  the 
interface  is  a  fraction  of,  and  opposite  in  sign,  to  that  at  the  top  of  the 
surface  layer,  in  the  first  law  of  thermodynamics,  this  latter  equation  is 
then  closed.  We  may  then  use  this  same  (well-known,  see  Carson,  1973) 
closure  in  the  interface  Equation  (14)  ,  together  with  the  various  modelling 

assumptions,  to  derive  a  prediction  equation  for  . 

oC 

The  ground  temperature  and  soil  moisture  equations  are  adapted  from 
Deardor ff  (1978). 

The  final  form  of  the  equations,  in  which  we  have  used  Leibniz's  rule 
in  the  form 


h  ^ d  *  r 


First  Equation  of  Motion 


-  4-  fl'  ^52  +  S~-  («  '-f  [M+i)^-/)^)n^ 

^  ^  rc*-0 

±  {[&[£)  _Atf+$)]^-  4[^4ftV«V?'- rtV-R+0]^^-  I 

.  .  o  *  ^  1^1  u 

U~^  “ 

Second  Equation  of  Motion 

2.  -+  ?rBar  +  ?3r  Stop¬ 
s'  11  izin 

xvUft')  +  [*ttL)Mk)jA({)'i{t)-'*t  -  A«-tsy«^)\s4- ' 

:  _  Q  >fQ  Ivl  V 

(*-V> 

First  Law  of  Thermodynamics 


W 


Inversion  Height  Equation 


I" 1  ^  -  6(i)  +£loVr]  (A --ttVr~  W]  * 

J  AU-K'l  &*!  -  Ai-fl^  -vft 

f<r>y  E«+S)&lt  -B«s)6*  +  ea-rt')'<^f«'>  i  «+S>W-««>'^ 


Dust  Concentration  Equation 


~  ^ \  p  ~^A  4  ^  u^i-  +  fcp  V^-  -  /u/f>-V-CL)l  r 

^«-4AL  -«  -w  *  i 


Q fe»  "Ca*) 
D  fJU-0 


Moisture  Equation 


^  ^*h«  ^<w,(u^U<  *  3n- 

*sMZtV£> 

t - (X40 - y*  <*•«•-  J  ~ 


4  M  g 


nl  *  cmo  vi 

5jS*t 

(A -A) 

Ground  Temperature  Equation 


>tr  6M.  t 

Ha  r  P4  Cho  A/-0  u  -  K were t 


Here  =  net  radiation  balance  at  the  ground 


Soil  Moisture  Equation 

=  -  c.  [IjjJI  .  <\ 

o  - -  ' 

C2  and  C2  are  constants  fw  P*  T, 

Pw  =  density  of  liquid  water 

W2  bulk  moisture  (analogous  to  O2) 

P  =  precipitation  rate  (prescribed  ' 

Ps  =  soil  density 

cg  =  soil  specific  heat 

dl  =  Ost!)172 

ks  =  soil  thermal  diffusivity 
^  =  period  of  1  day 

=  deep  soil  potential  temperature 

In  Equations  (50)  and  (53)  above,  there  appear  the  lapse  rates  of 

potential  temperature  and  mixing  ratio,  f(t),  and"?(t),  and  their 

derivatives  ^  and  ^  .  In  this  type  of  model,  it  is  not 

possible  to  derive  an  equation  for  prediction  of  these  quantities.  Thus  we 

have  specified  'f(t)  and$(t).  see  Fig.  2  for  the  form  of  )f(t). 


In  addition  to  the  prediction  scheme.  Equations  (48)-(55)  inclusive. 


there  exist  three  diagnostic  equations,  two  for  the  geostrophic  wind 
components,  and  the  equation  of  continuity  for  the  vertical  velocity. 


Geostrophic  Wind  Equations 

I;  W  +r 


Continuity  Equation 


The  expression  for  wk  is  derived  in  the  following  way: 
Substituting  Eq.  (20)  into  Eq.  (17) ,  we  have 


Then,  integrating  the  continuity  equation 


(60) 


I  --  (4+ 3.) A  fii^N -A 

~io  ' 

wT  is  given  by  Eq.  8,  and  in  that  equation  is  evaluated  from 


Given  initial  values  of  u,v,  9R,h,ck,qk,eGR,WGR,  the 
system  can  be  solved.  If  we  are  concerned  with  a  limited  area,  lateral 
boundary  conditions  are  important.  These  will  be  discussed  later.  The  upper 
and  lower  boundary  conditions  have  already  been  prescribed. 

If  we  assume  no  change  of  magnitude  or  direction  of  horizontal  wind  with 
height  within  the  transition  layer  only,  then  A(z)  =  B(z)  »  1.  This  situation 
freguently  exists  when  the  transition  layer  is  well  mixed.  The  equations  then 
become: 


First  Equation  of  Motion 


tt-1) 


(62) 


Second  Equation  of  Motion 


+  £-ra- 4  *  to- to + 

"rth  "W  '  0  77  »7 


+  fc-»««fWW>>£+ 


First  Law  of  Thermodynamics 


iM.  41(J  6^+^,  (v  6i)  +  tkA)  ht  4  Hk^ 


V  'Xt 


^  -A  —A.) 


+  -£  (mt«  +*0  -  4-  fn-ft.y^vKftta- 


Inversion  Heiqht  Equation 


u  IjU  «t4#l«+sy^(,>S  4«.A)^W»«Vr 

« ti)  6«  t  -  &4.  t  W  rt) W  Vjw*  «Vt  -  r]^  -fl ,  r»  \  v )  le^-af) 


Dust  Concentration  Equation 


IF  ^ +% (^)  +(^\^+sf 

=  Si  ill 

VU-A) 
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Moisture  Equation 


-V  -vfftr  +  ift)  U-4)  +  j£ 

^  9  ^  (*-&)  (*-4} 

+5u«IC:t^  -B«+iV-t-*M'i]v^;  +  ?^l1^  s]+^.  [c£J)-J-] 

-  e».  V*1  lv«  ,6^  ' 

The  Ground  Temperature  Equation  (54),  Soil' Moisture  Equation  (55),  and 
Geostrophic  Wind  Equations  (56)  are  unchanged. 

The  expressions  for  u  and  v  in  the  surface  layer  become 


u-- 


v  = 


Equations  of  Continuit 


u/;  _  ur 


From  Eq.  (69)  ,  we  deduce 


T 


-  K  ^)T’i 


This  expression  is  useful  in  a  model  in  which  w^  is  prescribed,  as  in 
the  one  -  dimensional  model  to  be  described  below.  For  a  range  of  values  of  h 
(300m*h*2000m) ,  k  =  20m,  ZQ  =  0.7cm,  and  with  w^  =  0,  we  find 

-.06  wh  t  wk  *-.01wh  ' 

r 

Thus,  in  this  model  w^  is  very  small,  and  is  essentially  zero,  unless 

=  0. 


Finally,  we  write  the  equations  for  a  model  in  which  horizontal  gradients 
are  absent. 


First  Equation  of  Motion 


<n> 


Second  Equation  of  Motion 


^  V©*  A)  TPu 


(72) 


First  Law  of  Thermodynamics 


iF  +0$  £(<*>**  VJi)  =  (r*'f£l  + 

- (T^~ 


Inversion  Height  Equation 


s  \AT^  _  cm>  (  8^  -  fyfe") _ 

(fefcl-  6^)+  %<•>-*]  ft-*) 


Dust  Concentration  Equation 


A-  iL.  -  (w.,  +-rvS\ .  Q  \v|(<jfc-cV>  (75) 
-xt  *  ;) 


Moisture  Equation 


—‘I*-  -y 

_3t  (4-4-') 

±fa-4.')\s-n.)\+S 

-fyu-feSi  3-t, 

*c~m*  CAL.  ( —  *?  aI  \kr 


Yv^ 
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The  Ground  Temperature  Equation  (54)  and  Soil  Moisture  Equation  (55)  are 


unchanged.  The  geostrophic  wind  Equations  (56)  cannot  be  applied. 


sou  ,v  must  be  specified.  Similarly,  the  Equations  of  Continuity 
'•-9  9 

(69)  cannot  be  applied,  so  that,  if  an  estimate  of  the  effect  of  w^  is 
desired,  w^  must  be  specified.  Then,  w^can  be*deduced  from  Eq.  (70). 


THE  FINITE  -  DIFFERENCE  SCHEME 


We  have  used  a  domain  300km  x  60km,  with4x  =  10km.  Fig.  3  shows  the 
Eliassen  grid  (Mesinger  and  Arakawa,  1976)  ,  which  is  a  space  -  time  grid 
staggered  in  both  space  and  time,  convenient  for  the  leapfrog  scheme  associated 
with  centered  space  differencing. 

Two  and  four  point  averages  of  variable  quantities  were  taken  as  needed  to 
provide  values  at  the  grid  point  under  consideration. 


Thus  averages  are  defined  either  as 


(77) 


or 


=  i 


(78) 


Centered  differences  were  taken  over  one  or  two  grid  intervals  as  dictated 


by  locations  of  variables  on  the  grid.  For  example,  when  predicting  u.^a 
term  like  (3h/&x)n  is  approximated  using  a  one-grid  interval  centered 
difference 


( ^  ’  A  ^  -j.  -  c-i-A 


(79) 


When  predicting  9.  .rtf^,a  term  like  Cbh/^x)  n  is  approximated 
1  rj 

using  a  two-grid  interval  centered  difference 


When  a  two  or  four  point  average  was  necessary  for  squared  terms,  the 
averaging  operator  was  performed  first,  then  the  average  was  squared. 

When  a  choice  was  necessary  between  one  or  two  grid  -  interval  differences 
for  nonlinear  terms,  such  as  (uO) ,  the  one  grid  interval  difference  was 
used.  All  products,  such  as  uO,  were  formed  after  appropriate  averaging  to 
provide  values  at  the  same  point. 


The  leap-frog  time  differencing  scheme  is 


LATERAL  BOUNDARY  CONDITIONS 


We  use  the  radiation  condition  proposed  by  Orlanski  (1976).  This  makes  use 
of  the  Sommerfeld  radiation  condition. 

r 


UL  +  C  =  O  (82) 

lit 7  ^ 


where  C  is  constant. 

The  finite  -  difference  leap-frog  representation  of  this  equation  is 


Here  JM  refers  to  the  boundary  point. 


The  essence  of  Orlanski' s  (1976)  technique  is:  instead  of  fixing  a 
constant  value  of  the  phase  velocity,  we  numerically  calculate  a  propagation 
velocity  from  the  neighboring  grid  points,  using  the  same  Equation  (83)  for 
each  variable  to  calculate  C.  Thus  we  find 


c  _  _  jV (JM-iV 


4>n  (j/i-I) 


L  <?"'  Pn-i) 


(84) 


Tb  determine  the  boundary  grid  point,  from  Equation  (84)  is  substituted  for 
C  in  Equation  (83)  ,  in  which  the  time  index  is  increased  by  1.  We  find 


In  this  formulation,  we  require  (XC^Ax/at. 


For  the  limiting  outflow  condition, 


=4X/4t, 


4mi  (5h)  =  4s  (srt-1) 

i 


(86) 


If  C^O, 

t*"  (Th)  =  4>""  (th) 


(87) 


No  information  has  come  from  the  interior  solution  when  =  0,  so  we  must 
regard  this  as  inflow  information,  to  be  prescribed  from  a  previous  time  step. 
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Finally,  we  have: 


Calculate  from  Equation  (84) . 


If  >  Ax/at  (limiting  outflow) ,  set  =  4x/at 
If  0  <  C^<  Ax/At  (outflow),  use  as  is 
If  <  0  (inflow),  set  C^=  0 


(S8> 


Use  from  Equation  (88)  to  calculate  <j>  (JM)  from  Equation  (85)  . 

In  our  initial  tests,  we  shall  use  only  the  condition  corresponding  to 
=  0,  i.e.,  Equation  (87),  on  all  boundaries,  since  this  is  an  easy 
condition  to  apply  for  testing  the  model  equations  themselves. 

In  discussing  the  above,  we  have  tacitly  assumed  that  we  are  dealing 

with  the  right  and  bottom  boundaries  respectively.  To  derive  the 

★ 

appropriate  equations  for  the  left  and  top  boundaries  ,  we  must 
We  have 


—  4^'Vlh) 


1 4"  -UM±±_ 

L  ^ 


(89) 


where,  due  to  the  method  of  indexing,  JW-1  now  refers  to  the  first  grid 
point  in  from  the  left  and  top  boundaries.  Applying  Equation  (89)  to 
neighboring  grid  points. 


* _ 

when  applied  to  the  top  and  bottom  boundaries,  a  x  is  replaced  by4y 


in  the  set  of  equations. 


Cf  -  -  L^Vxhti)  -  4>*"*  fTh+1)”] 

4VI  (~s h  (Th+i)  4-  4>h  * 


4*x 

O-A  t- 


(90) 


We  now  substitute  Equation  (90)  into  Equation  (89)  ,  using  for  C,  and 
increase  the  time  index  by  1. 


We  obtain 


In  this  formulation,  we  require  -Ax/At  <  <  0  . 

For  the  limiting  outflow  condition,  =  -Ax/At, 

=  4>"  (JH+t") 

'  (92) 

as  before. 

If  =  O  y 


4>*+'  (Trt)  =  d>*  ’  (tm) 


(93) 


This  is  again  regarded  as  inflow  information. 

Finally,  we  have: 

Calculate  from  Equation  (90) . 

If  <  -4  x/4t,  (limiting  outflow),  set  =  -  4x/dt 

If  -  4  xAt  <  0,  (outflow),  use  as  i&.  (94) 

If  >0  (inflow),  set  =  0 

Use  from  Equation  (94)  to  calculate  (jm)  from  Equation  (91). 

In  the  above  Equations  (84)  and  (90) ,  it  is  conceivable  that  may 
become  infinite  or  indeterminate. 


Case  1  -  Cx  infinite 


Equation  (85)  may  be  written 


(95) 


V-l  V*\ 


v\41 

+  (i  nN)= 


Similarly,  Equation  (91)  becomes 

Uw>  4  (in)  =  % 

C^-*>ioo 
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(97) 


Case  2  -  indeterminate 

This  situation  arises  when,  in  addition  to  the  fact  that  the 
denominator  of  Equation  (84)  is  zero,  the  numerator  is  also,  then 

<\n  (JM-1)  =4>n~2  WM-U  (98) 

If  applied  to  the  boundary  point  <JM) ,  this  is  essentially  the 

f 

statement  Equation  (87)  corresponding  to  =  0.  Thus,  in  case  the 
numerator  of  Equation  (84)  or  Equation  (90)  is  zero,  we  use  the  condition 
corresponding  to  =  0,  i.e..  Equation  (87)  or  Equation  (93). 


EXPERIMENTS 


One  -  Dimensional 


Our  initial  experiments  are  all  concerned  with  a  dry  atmosphere.  We 
first  ran  the  one  -  dimensional  model,  given  by  Equations  (71) -(75) 
inclusive  plus  Equation  (54) ,  with  the  following  conditions: 

A  (bt£)  =  B  (hf£)  =  1  (no  wind-shear  across  inversion) 

A  /* 

u  =  v  =  0  (no  geostrophic  wind) 

a.)  w^  =0;  =  0 


b.)  w,  =  -  (h-k) =  —  (h— k)  B  (B  =  constant);  w.  =0 


B  =  -5xlO"V( 0900-1 500) 

-5x10_5S-  +  5xl0”V('l500-2100) 
5x10“5j(2100-0300) 

5x10_5S-  -  5xl0-5*0300-0900) 


A,  A  —1 

c.)  u  *  v  «  100  cm  s 


d.)  ©k  =  303  °K,  eGR  =  303. 2°k 
••)  Cho  =  4  *  10"2 


£.) 

g. ) 

h. ) 

i. ) 


h  =  300  m 
o 

Ck  =  100  /<gnf  3 


C_  =  100  jugm  =  constant; 

CjR 

C  =  10,000  /wgm-3  =  constant 


The  results  of  the  experiment  with  wh  =  0  are  shown  in  Figures  4,  5, 
and  6. 

A  A 

Figure  4  shows  the  evolution  of  the  fields  of  u  and  v  (the  caps  are 
omitted  from  the  figure) .  These  fields  evolve  smoothly.  The  absolute  value  of 
the  wind  reaches  a  maximum  at  about  1300  LST  and  a  minimun  at  about  0100  LST. 
Fig.  5  shows  the  evolution  of  the  ground  temperature  On* and  the  temperature  at 
the  top  of  the  surface  layer  @k.  0^  reaches  its  maximum  value  at  about 
1400,  while  0k  reaches  its  maximum  at  about  1600  local  time.  Both  of  these 
are  reasonable.  The  ground  temperature  reaches  its  minimum  at  about  0600, 
while  the  air  temperature  reaches  its  minimum  at  about  0800,  which  is  somewhat 
late.  The  maximum  and  minimun  of  both  exceed  the  corresponding  values  for, 
which  is  realistic.  The  range  of  both  values  is  somewhat  less  than 
should  be  expected  on  a  summer  day  with  light  winds. 

Figure  6  shows  the  evolution  of  the  inversion  height  (H  in  the  figure)  and 
dust  concentration  for  the  two  cases. 
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CGR  =  lOOyUgm-3,  CQR  =  10,000/<.gnf3. 

The  former  of  these  corresponds  to  an  average  air  value  for  Saharan  dust.  It 
turned  out,  after  6  months  of  observation  at  Sde  Boqer,  that  the  average  value 
is  closer  to  60/jgm-3.  However,  the  use  of  100  will  not  alter  these 
results  in  any  substantial  way.  The  value  10,000/Agm-3  corresponds  to  that 
within  a  strong  dust  storm.  The  first  thing  that  can  be  seen  from  this  figure 
is  that  the  inversion  height  rises  steadily  throughout  the  day.  This  may  be 

f 

indicative  of  the  fact  that  a  convective  boundary  layer  model  should  apply  only 
during  daylight  hours  -  or  it  may  indicate  that  subsidence  at  inversion  level 
is  required  to  bring  the  inversion  down  at  night. 

The  curve  CGR  =  100  yM  gm- 3  shows  the  trend  of  (dust 
concentration  at  20  meters).  This  quantity  decreases  steadily  with  time.  This 
is  not  surprising,  since,  according  to  Equation  (75),  the  tendency  of  is 
inversely  proportional  to  the  inversion  height.  It  also  indicates  that  the 
turbulent  transport  of  dust  was  too  weak  to  overcome  the  thinning  out  of 
as  the  inversion  rose.  With  CGR  =  10,000yugm~3 ,  there  is  continuous 
intense  turbulent  transport,  so  that  the  dust  concentration  at  z  =  k  is  kept 
high  throughout  the  period. 

It  should  be  recalled  that  both  ek  and  h  are  dependent  upon  the  assumed 
form  of  'f(t)  (Figure  2).  Thus,  it  is  possible  to  "tune"  the  model  by 
experimenting  with  other  curves  of  T(t) . 

Figures  7,8,9  shows  the  results  for  the  same  quantities  when  w^,  the 
vertical  velocity  at  inversion  height,  is  not  zero.  In  actual  fact,  this 
quantity  should  be  obtained  by  solving  the  continuity  equation,  which  is  not 
available  in  a  model  without  horizontal  transport.  Thus  we  assume  values  of 
the  horizontal  divergence  B  (given  in  the  list  of  data) ,  which  vary  with  time 
throughout  the  day.  Again,  this  quantity  is  "tunable",  but  we  wish  to 
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highlight  the  effect  of  subsidence  on  inversion  height. 

Figure  7  shows  the  evolution  of  u  and  v  when  w^  ^  0.  The  results  are 
indistinguishable  from  those  of  Figure  4.  This  is  simply  because  the  w^ 

terms  vanish  from  Equations  (71)  and  (72)  under  the  conditions  we  have  assumed. 

Figure  8  shows  the  evolution  of  ^GR  and  0^.  There  are  some  small 
differences  between  these  results  and  those  with  w^  =  0,  but  the  evolution 
is  more  or  less  the  same  in  both  cases.  r 

The  major  differences  in  the  case  wh  f  0  can  be  seen  in  Figure  9,  when 
compared  with  Figure  6.  Figure  9  shows  the  evolution  of  inversion  height,  h, 
and  of  dust  concentration,  C^.  The  inversion  height  now  reaches  a  maximum 
of  about  1400  m  at  about  1800  LST,  and  then  falls  due  to  subsidence.  It  reaches 
a  minimum  of  about  600  m  at  about  0500,  and  then  starts  up  again.  Even  though 
the  value  itself  at  the  minimum  may  not  be  realistic,  this  experiment  very 
strongly  highlights  the  effect  of  the  vertical  velocity  at  inversion  height 
level  on  the  inversion  height  itself. 

The  same  type  of  result  is  true  for  C^.  For  CGR  *  100/Kgm-^, 
there  is  little  change  in  C^,  except  that  it  does  not  decrease  as  much  as 
it  did  in  the  case  w^  =  0.  This  is  so  because  the  subsidence  during  the 
latter  part  of  the  period  tends  to  increase  the  concentration  at  lower  levels. 

This  effect  is  even  more  marked  in  the  case  of  CGR  =  10,000/Ugm-3. 

There  is  intense  turbulent  transport  upward  during  the  entire  period.  During 
the  second  half  of  the  day,  the  upward  transport  to  level  z=k  is  intensified  by 
subsidence  to  that  level,  so  that  continues  to  increase. 

It  should  be  realized  that  the  assumption  of  CGR  =  constant  for  24 
hours  is  rather  weak,  in  the  sense  that  wind  erosion,  which  causes  variations 
in  CGR,  occurs  in  bursts.  We  shall  discuss  this  point  later. 

The  legends  at  the  top  of  all  the  figures  contain  the  heading  "NCN  - 
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PARAMETERIZED  DUST  EQ."  We  have  distinguished  these  results  from  results  with 


"PARAMETERIZED  DUST  EQ" ,  to  be  shown  below,  for  the  following  reasons. 

We  have  used  the  interface  Equations  (12)— (16)  inclusive  ,  to  eliminate  the 
turbulent  fluxes,  at  z=h  except  for  Equation  (14) .  It  is  also  possible  to 
parameterize  the  other  fluxes,  for  example  (w'c')h  in  Equation  (16).  If 
we  use  the  flux  -  gradient  relation 


-w'c' 


=  A*(z) 


(99) 


then  -(w'c')h  =  A*(h)  D'(h)ck 


(100) 


where  A  (z) ,  the  coefficient  of 


eddy  diffusion  ,  is  given  by  Equation  (46) . 


In  that  case.  Equation  (16)  becomes 


^ fr 


Df-C) 


(101) 


and  the  dust  concentration  equation  becomes 


-  OVvKc^-cV)  f  ft*  Ql 
l  (X  -4) 

instead  of  Equation  (75) .  We  now  see  that  the  model  in  this  form  contains  two 


inversion  height  equations,  Equation  (74)  due  to  convection  of  heat  and 
Equation  (102)  due  to  turbulent  diffusion  of  dust.  We  carried  out  experiments 
with  the  same  initial  data  ,  but  in  which  the  inversion  height  was  calculated 
as  an  average  of  the  results  from  Equation  (74)  and  (102) .  The  results  are 
shown  in  Figures  10-15  inclusive.  The  major  differences  between  these  results 
and  the  non-parameterized  results  are  in  the  inversion  heght  and  dust 
calculations,  seen  in  Figures  12  and  15.  These  should  be  compared  with  Figures 

6  and  9  respectively.  For  the  w^  =  0  case,  the  results  show  similar  trends, 

but  the  non  -  parameterized  inversion  height  went  much  higher,  and  the  curve  is 

less  smooth.  This  is  undoubtedly  due  to  the  averaging  involved  in  obtaining 

the  parameterized  height.  The  variation  of  is  quite  similar  in  both 
cases,  although  when  CGR  =  10,000 jUgm  does  not  drop  down,  in 

the  parameterized  case  as  it  did  in  the  non-parameterized  case.  When  w^  f 
0  the  curves  look  quite  similar,  except  that  the  maxima  and  minima  of  inversion 
height  occur  about  an  hour  later  in  the  parameterized  case  than  in  the 
non-parameterized  case. 

The  conclusion  to  be  drawn  from  these  comparisons  is  that  the 
non-parameterized  system,  which  is  theoretically  more  correct,  should  be  used. 


Two  -  Dimensional  Experiments 

We  used  the  system  of  Equations  (62)  -  (66)  inclusive.  Equations  (68) , 
(69),  (70).  We  have  not  assuned  A  (htj),  B  (ht$)  =  0,i.e.,  we  have  allowed 
wind  shear  through  the  inversion. 


a) .  Constant  Boundary  Conditions 

We  consider  a  region  300km  x  600km,  with  a  grid  spacing  4  x  =  dy  =  10km. 
this  experiment,  we  used  the  lateral  boundary  condition  corresponding  to  C^ 
0,  i.e., 

(j.Y\)  =  4>"~'  fori)  187) 

i 

r 

Thus  we  expect  that  this  fixed  boundary  condition  will  eventually  lead  to 
instability  due  to  reflection  at  the  boundaries. 

As  this  experiment  was  essentially  a  check  of  the  program,  we  used  very 
simple  initial  conditions. 

/v  —1 

u  =  100  cm  s  everywhere 

A 

v  =0 

h  =  650  m 

6k  =  303°k 

6gr  =  303. °2k 

)f(t)  as  in  Figure  2. 

Ck  =  100/Ugm  -3 
CGR  =  100/lgm-3 

A(ht$)  =  B(h+S)  =  1.01  -  corresponding  to  1ms-3  (10m)-3 

A  _1 

Ug  was  taken  as  1  ms 


In 


We  expect  changes  in  the  variables  due  to  changes  in  the  radiation  fluxes. 


We  used  a  time  step  of^t  =  1  minute.  When  the  proqram  ran  smoothly 
for  33  time  steps,  we  tried  4t  =  2  minutes,  with  identical  results.  We 
then  ran  for  a  lonqer  period,  but  the  calculations  blew  up  at  time  step  122 
*  4.07  hours,  primarily  due  to  boundary  instability.  The  evolution  of  the 
interior  fields  seemed  reasonable.- 

f 

b. l  Radiation  Boundary  Conditions 

Startinq  from  the  same  initial  fields,  we  applied  the  radiation 
conditions  Equations  (821  -  (981  inclusive.  The  experiments  have  been  run 
®o  far  for  6  hours,  with  one  minute  time  steps.  The  evolution  of  the 
interior  fields  seemed  reasonable.  Stronq  qradients  did  not  develop  in  all 
f-he  interiors  as  they  did  in  the  constant  boundary  condition  case. 

c. l  Comparison  of  Results  of  the  Two  Experiments 

In  the  fiqures  to  be  discussed  below,  we  have  displayed  the  printout® 
for  the  various  fields,  only  for  rows  1,2,3.4,5.58,59,60,61.  The  remaininq 
fields  are  smooth  transitions.  Due  to  the  printinq  limitations,  each  row 
of  the  qrid  is  represented  by  almost  two  full  rows  on  the  printout,  so  that 
columns  1  throuqh  17  qo  from  left  to  riqht,  while  columns  18  throuqh  31, 
just  beneath  them,  qo  from  riqht  to  left.  We  have  since  modified  the  print 
routine  to  qive  a  clear  rectanqular  array,  and  future  results  will  b« 
printed  in  the  manner. 

Piqure  16  shows  the  u-field  at  180  minutes,  with  constant  boundary 
conditions.  Fiq.  17  shows  the  Vfield  with  radiation  boundary  conditions. 
Althouqh  the  nunbers  on  the  bottom  boundary  of  the  latter  are  larqer  than 
in  the  constant  B.C.  case,  the  fields  everywhere  else  are  quite  uniform  and 


smooth.  Fiqures  ]8  and  19  show  the  v-fields  in  the  constant  and  radiation 
cases.  Aqain,  the  interior  fields  in  the  radiation  case  are  much  smoother, 
and  less  subiect  to  qradients  created  by  boundary  reflection.  The 
inversion  heiqhts  are  shown  in  Fiqures  20  (constant)  and  21  (radiation) . 
There  is  no  question  of  the  superiority  of  the  radiation  condition  for  this 
field,  even  thouqh  the  left  and  upper  boundaries  are  not  sufficiently 
smooth.  Fiqures  22  and  23  show  the  fields* in  the  constant  and 
radiation  cases.  Except  for  the  bottom  row.  the  radiation  condition  qives 
a  smoother  field.  Fiqures  24  and  25  show  in  the  constant  and 
radiation  cases.  Aqain.  the  radiation  case  is  clearly  superior.  Fioures 
•56  and  27  show  the  (dust)  fields  in  the  constant  and  radiation 
cases.  In  this  field,  the  interior  for  the  constant  case  alredv  shows 
larqe,  unrealistic  values,  while  these  values  for  the  radiation  case  are 
very  smooth.  Finally,  Fiqures  29  and  29  show  the  wh  values  and 

A  A 

radiation  cases.  These  values  should  all  be  zero,  as  both  u  and  v  should 
be  developinq  uniformly.  While  all  are  very  small,  the  fields  in  the 
constant  case  are  already  beqinninq  to  develop,  while  they  are  almost  all 
zero  in  the  radiation  case. 

These  results  clearly  illustrate  the  efficacy  and  superiority  of  the 
radiation  boundary  condition  over  the  constant  boundary  condition.  Yet  the 
boundaries  still  require  additional  treatment.  We  plan  to  experiment  with 
suitable  filters  (Shapiro.  1975)  to  control  spurious  hiqh  frequency 
oscillations  on  the  boundaries  themselves.  There  are  also  suqqestions 
(Miller  and  Thorne,  1981)  for  improved  versions  of  the  radiation  boundarv 
condition. 

As  noted  earlier,  the  radiation  boundarv  condition  experiments  have 
been  continued  for  6  hours,  with  one  minute  time  steps.  All  of  the  fields 
developed  reasonably,  with  little  encroachment  from  the  boundaries,  except 
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the  inversion  heiqht  field.  While  the  calculations  still  did  not  blow  ud 
at  6  hours,  spurious  oscillations  were  developing  in  this  field.  It  is 
clear  that  some  kind  of  smooth inq  is  required.  This  will  be  attempted  in 
future  experiments. 


CONCLUSIONS 


Ttie  one  -  and  two  -  dimensional  versions  of  the  desert  planetary 
boundary  layer  model  appear  to  be  checked  out,  in  the  sense  that  they  qive 
reasonable  metereoloqical  results. 

The  one  -  dimensional  version  runs  well  for  24-hours,  and  can  certainly 
be  run  for  lonqer  periods.  It  is  now  possible  to  use  this  model  for  a 
rvunber  of  sensitivity  experiments,  such  as  prescription  of  dust 
concentration  at  the  qround  as  a  function  of  time,  variable  vertical 
velocity  at  inversion  heiqht,  variable  surface  albedo. 

The  two  -  dimensional  version  runs  well  for  up  to  6  hours  before 
fluctuations  in  the  inversion  heiqht  field  set  in.  Thus,  more  nimerical 
investigation  is  needed  in  order  to  carry  the  experiments  further.  As 
already  stated  above,  the  investigations  will  take  the  form  of  smoothinq 
and  filtering  operators,  and  possible  improved  forms  of  the  radiation 
boundary  conditions. 

We  have  already  prepared  a  set  of  base  data  of  all  fields  for  Israel, 
and  will  carry  out  experiments  with  these  data,  over  irreqular  terrain, 
following  adoption  of  more  suitable  numerical  techniques. 
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RECOMMENDATIONS 


There  are  a  number  of  steps  which  can  be  taken  which  may  lead  to 
improvement  of  the  present  model.  Although  the  model  appears  to  be  checked 
out,  in  the  sense  that  it  does  not  blow  up  by  6  hours,  and  gives  reasonable 
evolution  of  the  fields  during  that  time  it  is  desirable  to  focus  on  a 
model  which  may  be  useful  operationally. 

There  are  still  several  steps  which  may  be  taken  to  improve  the 
radiation  boundary  conditions.  Among  than  are  the  application  of  a 
suitable  filter  (Shapiro, 197 5) ,  which  was  successfully  applied  to  this  type 
of  problem  by  Eliassen  and  Thorsteinsson  (1984).  Another  is  the 
possibility  of  applying  an  improved  version  of  the  radiation  boundary 
condition  according  to  the  method  suggested  by  Miller  and  Thorpe  (1981). 

There  are  also  several  ways  of  improving  the  physics  in  the  dry 
model.  We  could  try  a  possibly  more  realistic  version  of  the  stability 
Y(t) .  We  could  also  try  to  improve  the  assumption  that  the  dust 
concentration  near  the  ground,  C^,  is  constant.  To  do  so  requires  a 
method  for  predicting  erosion  of  soil  by  wind  .  A  first  step  in  this 
direction  has  been  taken  by  Berkofsky  (1984)  who  developed  an  equation 
describing  the  processes  of  detachment,  transport  and  deposition  which  make 
up  wind  erosion. 

When  all  of  the  above  has  been  incorporated  we  can  then  include 
moisture  and  a  thermally  active  surface  layer.  Real  data  (already  compiled 
for  Israel  for  this  model)  can  then  be  used  as  input,  together  with 
topography,  for  operational  testing. 
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DUST  CONCENTRATION  DATA 


In  Table  1  we  present  dust  concentration  data  for  Sede  Boqer,  compiled 
simce  June  1984.  For  a  variety  of  technical  reasons,  the  data  are  not 
continuous.  These  dust  concentration  data  were  obtained  with  an  ultra  - 
high  volume  sampler  (Sierra  Instruments)  with  cascade  impactor,  with  a 

f 

constant  flow  meter  operating  at  40  cubic  feet  per  minute.  These 
categories  are: 


C1 

7.2-00  /M. 

C2 

3.0  -  7.2  M 

C3 

1.5  -  3.0  M 

C4 

0.95  -  1.5  yU 

C5 

0.49  -  0.95 

C6 

<  0.49  yU 

_3 

The  averages  in  ^igm  ,  are: 

Total  C:  C2  C3  C4  C5  Cg 

59.81  11.24  15.23  12.13  10.26  5.98  18.71 


The  minimum  value  was  10.47 ^ignf^  on  11/11/84.  The  maximum  value 
was  133.  20/<gm  on  1/8/84.  The  wind  direction  is  given  as  a  sort  of 
twenty  -  four  hour  average.  It  is  seen  that  westerly  winds  dominate. 

These  data  will  have  to  be  analyzed  with  respect  to  the  prevailing 
synoptic  situation. 
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Table  1.  -  IXist  concentration  Data  at  Sede 
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Fig.  19.  -  v  field  at  180  minutes,  radiation  B.C 
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Fig.  21.  -  Inversion  height  at  180  minutes,  radiation  B.C. 
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Fig.  27.  -  at  180  minutes,  radiation  B.C. 
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